function [x,rn,rR,reps] = romberg(f, a, b, eps, nmin, nmax, mmax, debug)
% Aufruf x = romberg(f, a, b, eps, nmin, nmax, mmax, debug) 
% f ... Funktion
% a ... Intervallanfang
% b ... Intervallende
% eps ... realtive Fehlertoleranz
% nmin ... 2^nmin minimale Anzahl der Auswertungspunkte
% nmax ... 2^nmax maximale Anzahl der Auswertungspunkte
% mmax ... maximale Anzahl der Extrapolationsschritte
% debug ... Wenn debug auf true gesetzt ist, gibt die Funktion romberg fuer Experimentierzwecke viele Informationen ueber den Verlauf der Approximation zurueck.
% x ... Approximation des Integrals

% Da unter anderem auch die Qualitaet der Approximation und die Hintergruende fuer Qualitaetsunterschiede untersucht werden sollen, kann romberg durch Uebergabe des debug-Parameters dazu aufgefordert werden viele weitere Informationen auszugeben. Unter anderem laesst sich auch die Speicherverwaltung so umstellen, dass alle Approximationswerte und nicht nur die letzten gespeichert werden.
% All dies ist zur normalen Verwendung zur Integration von Funktionen nicht notwendig und die Abfragen ob debug aktiviert ist, koennen sogar die Leistung negativ beeinflussen. Daher sind die Codestellen, die nur zu Debugzwecken dienen und entfernt werden koennen, wenn Debugging nicht gewuenscht ist durch die folgenden Marker markiert. (In C waeren das macros, aber das gibt's in Matlab nicht.)
%<--

%-->

rho=1; %Sicherheit eps

%<--
if (nargin < 8), debug=false; end
%-->

if (nargin < 7), mmax=nmax; end


%<--
global DEBUGROMBERG
DEBUGROMBERG=debug;
if(DEBUGROMBERG==true)
    R=zeros(nmax+1,nmax+1);
else
%-->

    R=zeros(2,nmax+1);

%<--
end
%-->

h=b-a;
R(ind(1),1) = h*(f(a)+f(b))/2;

for n=2:nmax+1
    h=h/2;
    R(ind(n),1) = R(ind(n-1),1)/2+h*sum(map(f,(a+(2*(1:2^(n-2))-1)*h))); %Wenn die Funktion f (wie die meisten Standardfunktionen in Matlab) sowieso vektorisiert ist, dann kann hier auf die Verwendung von map verzichtet werden.
    mgr=min(n,mmax+1);
    for m=2:mgr;
        R(ind(n),m)=R(ind(n),m-1)+(R(ind(n),m-1)-R(ind(n-1),m-1))/(4^(m-1)-1);
    end
    if n>nmin && abs(R(ind(n),mgr)-R(ind(n-1),min(n-1,mmax+1)))<rho*eps
        break;
    end
end

x=R(ind(n),mgr);

%<--
if(DEBUGROMBERG==true)
if (nargout >= 2), rn=n-1;
if (nargout >= 3), rR=R;
if (nargout >= 4), reps=R(ind(n),mgr)-R(ind(n-1),min(n-1,mmax+1)); 
end;
end;
end;
end;
%-->

end

function i = ind(n)
%<--
global DEBUGROMBERG
if(DEBUGROMBERG==true)
    i=n;
else
%-->

    i=2-mod(n,2);

%<--    
end
%-->

end

